data = readtable('../../data/main_VAR/UK_source_pre1946.xlsx', 'Sheet', 'Full_Data_for_VAR');
date = data.Year;

startdatenum = 1729;
enddatenum = 1946;

startdate = find(date == startdatenum);
enddate = find(date == enddatenum);
T_post = length(date(startdate:enddate));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

option.figure = 0;

Tstart = startdate; 
Tend = enddate; 
T = Tend - Tstart + 1;

date = date(startdate:enddate);

data = data((data.Year >= startdatenum - 1), :);
data = data((data.Year <= enddatenum), :);

taxrevgdp = data.revtogdp(2:end); % government tax revenue to GDP ratio
taxrevgdp0 = data.revtogdp(1); % government tax revenue to GDP ratio
spendgdp = data.spendingtogdp(2:end); % government spending before interest exdp. to GDP ratio
spendgdp0 = data.spendingtogdp(1); % government spending before interest exp. to GDP ratio
surplusgdp = taxrevgdp - spendgdp; % primary surplus to gdp
rpcgdpgr = data.realGdpGrowth(2:end); % real gdp growth (in logs)
inflation = data.inflation(2:end); % growth in GDP price deflator (in logs)
ynom1 = data.short_r(2:end) / 100; % nominal yield on a 1-year Treasury bill, expressed per annum
ynom10 = data.rate_10year(2:end) / 100; % nominal yield on a 10-year Treasury bond, expressed per annum

ynom1 = log(ynom1 + 1);
ynom10 = log(ynom10 + 1);

dy = data.dy(2:end);
pdm = log(1 ./ dy); % log price/dividend ratio
divgrm = diff(log(data.dy) + data.stockindex); % log dividend growth

deltalogg = [log(spendgdp(1)) - log(spendgdp0); log(spendgdp(2:end)) - log(spendgdp(1:end - 1))];
deltalogtau = [log(taxrevgdp(1)) - log(taxrevgdp0); log(taxrevgdp(2:end)) - log(taxrevgdp(1:end - 1))];

yieldspr = ynom10 - ynom1;
y0nom_1 = mean(ynom1);
y0nom_10 = mean(ynom10);
yspr0 = mean(yieldspr);
pi0 = mean(inflation);
x0 = mean(rpcgdpgr);
tau0 = mean(deltalogtau);
g0 = mean(deltalogg);
surplus0 = mean(surplusgdp);

% growth of log div to gdp ratio
divgrm = divgrm - inflation;

A0m = nanmean(pdm);
% mu_m = nanmean(divgrm);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Detrend tax and spending
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



tts(1) = taxrevgdp(1);
gts(1) = spendgdp(1);

for t = 2:T
    gts(t) = gts(t - 1) * exp(deltalogg(t) - g0);
    tts(t) = tts(t - 1) * exp(deltalogtau(t) - tau0);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Read in debt/gdp
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
gdebt = data.debttogdplevel;
logdts = log(gdebt(2:end));

deltalogd = [0, (logdts(2:end) - logdts(1:(end - 1)))'];

deltalogd_demean = deltalogd - mean(deltalogd);

run ../../tools/run_var_estimate_debtinvar;

save(['MAT\res_annual_VAR_debt_4CI_UK.mat'],'-regexp', '^(?!(option|Globaloption)$).');
